Liquid, Glass and Crystal in Two-dimensional Hard disks 
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We study the thermodynamic and dynamic phase transitions in two-dimensional polydisperse 
hard disks using Monte Carlo methods. A conventional local Monte Carlo algorithm allows us to 
observe a dynamic liquid-glass transition at a density p G , which depends very little on the degree of 
polydispersity. We furthermore apply Monte Carlo methods which sample the Boltzmann equilib- 
rium distribution at any value of the density and polydispersity, and remain ergodic even far within 
the glass. We find that the dynamical transition at p is not accompanied by a thermodynamic 
transition in this two-dimensional system so that the glass is thermodynamically identical to the 
liquid. Moreover, we scrutinize the polydispersity-driven transition from the crystal into the disor- 
dered phase (liquid or glass). Our results indicate the presence of a continuous (Kosterlitz-Thouless 
type) transition upon increase of the polydispersity. 



For several decades now, enormous interest has been 
concentrated on the glass transition in classical liquids 
|J. At this transition, the viscosity of a structurally dis- 
ordered liquid increases dramatically and reaches values 
typical of a crystalline solid. The underlying reasons for 
the spectacular slowdown of the system's time evolution 
have been hotly debated and views very much diverge on 
this question |g-^| . The transition is naturally related to 
the loss of thermal equilibrium during the quench towards 
low temperatures, and it is well-known that the proper- 
ties of a glass strongly depend on the detailed quench 
procedure |(J. The theoretical description of this phe- 
nomenon in all its complexity is a very difficult task. 

Notwithstanding the intimate relationship of the glass 
transition with loss of thermal equilibrium (i.e. ergodic- 
ity breaking), it is nevertheless possible to describe the 
glass within thermal equilibrium. This means that one 
studies configurations of the system taken from the Boltz- 
mann distribution (which exists in a mathematical sense 
at any temperature and pressure) and computes their 
thermodynamic and dynamic properties. For any ther- 
modynamically stable glass, the reference to the quench 
history can thus be omitted, and the theoretical analysis 
is much simpler. 

The most serious problem with the 'equilibrium glass' 
is of course that the Boltzmann distribution is experi- 
mentally inaccessible inside the glassy phase. Attempts 
to characterize the true equilibrium properties of a glass 
necessarily rely on the extrapolation of experimental 
data. 

However, computational methods have been developed 
which allow to bypass this obstacle: non-local Monte 
Carlo (MC) algorithms can remain ergodic in the glassy 
phase for prominent model systems, such as hard spheres, 
and allow to compute their equilibrium observables di- 
rectly pj . With this approach, we may address questions 



at the heart of long-standing controversies, for example 
about the existence of an inaccessible thermodynamic 
phase transition within the glass, which has been ren- 
dered responsible for the dynamic transition |2|-|5| . 

In the past, many classical statistical physics models 
have been considered in the study of glasses. For a long 
time, it was believed that the very simplest model sys- 
tems (monodisperse hard spheres, hard disks, and cen- 
tral potential models) did also possess a glassy phase 
j8[. However, more recent MC simulations have shown 
that monodisperse hard spheres always crystallize, and 
that the previously observed behavior was mainly due 
to finite-size effects and limited computer resources M. 
This has also been confirmed in experiments with col- 
loidal systems, where crystallization preempts the glass 
transition for monodisperse hard spheres. In order to 
avoid parasitic effects, some of the experiments have been 
performed in zero-gravity environments flO] . 

The situation changes abruptly, both in experiment 
flllfl and in the simulations, in the presence of a small 
dispersity in size. Several investigations have shown the 
existence of a so-called 'terminal polydispersity' of a few 
percent (c/ ; e.g., ||12| , |l3| ), above which crystallization can 
be avoided for any value of the external pressure, i.e., up 
to very high densities. 

Polydispersity opens up opportunities to study two 
very interesting, yet distinct phenomena: (i) the glass 
transition, which is no longer preempted by crystalliza- 
tion at sizeable values of the polydispersity; (ii) the ther- 
modynamic transition of the crystal into a disordered 
state (liquid or glass) upon an increase of the polydisper- 
sity. 

In the present work, we study these two phenomena in 
two-dimensional hard disks. We investigate their dynam- 
ical and thermodynamical phase diagram as a function of 
density p (or pressure P) and of the polydispersity for a 
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given functional form of the size dispersion. At polydis- 
persities for which crystallization is impossible, we find 
a purely dynamic liquid-glass transition at a density p , 
which depends very little on the degree of polydispersity 
e. For moderate values of e, the glass transition density 
can be defined even better than in the highly polydis- 
perse case considered earlier M . Interestingly, p G is con- 
siderably higher than the crystallization density for the 
monodisperse system. 

Our vastly improved equilibrium simulation methods 
allow us to scan the complete parameter space (e, p). 
Even far inside the glassy phase, we can compute the 
equation of state and the compressibility with very high 
precision. The polydispersity-driven transition of the 
crystal into the disordered phase (liquid or glass) also 
turns out to be particularly interesting: we find strong 
evidence that this disorder-driven phase transition is con- 
tinuous. Our findings can hardly be reconciled with the 
presence of a first-order transition, and differ from what 
was stated in a closely related work Q . 

In principle, for a polydisperse mixture of particles, one 
should choose the particle sizes ri independently from a 
given probability distribution p(r) . We would then have 
to perform an ensemble average over the distribution 
of radii. In dynamical simulations, this average corre- 
sponds to ensemble-averaging time correlation functions 
of a given sample. In thermodynamic calculations, the 
average over radii can be incorporated directly into the 
MC evaluation of the partition function, as some authors 
have done [Q. 

We have adapted a much simpler, yet practically equiv- 
alent approach, by considering a 'fixed probability incre- 
ment' model for any probability distribution p(r) with 
r > 0, where we pick the N radii rj e.g. according to the 
following equation 



drp{r) 



N+l 



(1) 



(with L drp{r) = 1). In this way, sample-to-sample 
fluctuations are completely eliminated, and the most rep- 
resentative sample is generated for an arbitrary distribu- 
tion pir). 

Here, we study a model with a flat distribution of radii: 



p{r) 



const r r , 
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otherwise, 



(2) 



i.e. the radii are given by r^ = r% + A(i — 1) (with 
A = (r„ - ri)/(N - 1)). As in Q, we set the total par- 
ticle volume as V par t = Nit /A and the Boltzmann factor 
as j3 — 1. 

We define the polydispersity by the normalized width 
of the distribution eq.(0) 



Notice that the distribution eq.(||) satisfies e < I/a/3, 
this limit being attained for r m i n = 0. 

Our simulations are performed in a periodic box of size 
L x \/3L/2, as is customary pa], and the only difference 
in set-up between our thermodynamic and dynamic sim- 
ulations is the choice of ensembles. The equation of state 
V(P) is most naturally computed in the [NP] ensemble, 
where the size of the box is allowed to fluctuate. The 
dynamical calculations are done in the [NV] ensemble, 
where the time evolution of states can be best monitored. 
For these calculations, a conventional local MC algorithm 
is used. Up to a trivial rescaling of time, we expect this 
algorithm in the long-time limit to give the same results 
as a molecular dynamics method. 

The equilibration of polydisperse systems is a difficult 
task, because the properties of the system depend sen- 
sitively on the relative arrangements of particles, which 
have to be averaged over. Therefore, one needs an algo- 
rithm that moves particles over large distances and av- 
erages over different possible arrangements. In the fluid 
phase this can in principle be achieved by a long series 
of local moves. In the crystal or the glass, it is however 
mandatory to equilibrate the system by non-local moves, 
as in the pivot cluster algorithm p6[ . 

A polydisperse system with a continuous distribution 
of radii is easier to simulate than a binary distribution, 
because it is often possible to exchange particles with 
similar radii. These swaps also implement non-trivial 
long-range moves, and have to be combined with a lo- 
cal Monte Carlo algorithm. The acceptance probability 
for such swaps decreases only at the most extreme densi- 
ties, when each particle is so close to its neighbors that it 
can virtually never be replaced by a particle with slightly 
larger radius. 

In high-density binary mixtures (c/ [|17|), the pivot 
cluster algorithm is clearly more appropriate, because 
the direct particle exchanges freeze out. In the present 
polydisperse case, however, we have obtained identical re- 
sults with both algorithms for systems of N = 256 and 
N = 1024 particles even slightly below what we believe 
to be close packing. 

Extensive tests of ergodicity (cf M) fully confirm the 
above statements: In very small systems (iV = 15), where 
the minimal size differences between particles are larger, 
we can show that the direct swap algorithm falls out of 
equilibrium at slightly lower densities than the pivot clus- 
ter method for the above-mentioned reasons. 

We point out a completely unrelated, yet crucial, issue: 
To compute the equation of state at constant pressure, we 
need to include volume changes into the MC algorithm. 
Extremely time-consuming sampling of the volume by 
trial and error (using the Metropolis algorithm |ll|) can 
be avoided 
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e = \J< Sr 2 



> J < r > 2 



(3) 



For orientation, we show in figure m the schematic 
phase diagram as a function of density p and polydisper- 
sity e (cfeq.(pl)) which results from our thermodynamic 



and dynamic calculation, and which we discuss in the 
remainder of this paper. We have studied this diagram 
by the routes indicated by gray lines in the figure, using 
either the non-local MC algorithm (for thermodynamic 
averages), or the local MC method (to study the dynam- 
ics). 
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FIG. 1. Schematic phase diagram of the two-dimensional 
hard-disk system as function of polydispersity e and density p. 
The dark lines indicate estimated phase boundary (obtained 
by linear interpolation of data points). The gray lines corre- 
spond to the paths along which calculations were performed. 
Full lines: thermodynamic calculations (non-local MC algo- 
rithm); broken lines: dynamic calculations (local MC). 

In Figure |2ja, we show the liquid and solid branches 
of the thermodynamic equation of state as a function of 
pressure P for a fixed e = 0.052. The presence of a ther- 
modynamic transition is evident. 
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FIG. 2. Inverse compressibility k _1 as a function of den- 
sity p and equation of state V(P) for two values of the poly- 
dispersity e at N — 256. Left: e — 0.052: a thermodynamic 
phase transition is clearly detectable, both in V(P) and in 
k" 1 . Right: e = 0.19: No thermodynamic phase transition 
can be detected, k is obtained by derivation of the curve 
V(P) and by direct computation from the fluctuations of the 
volume. 



In these calculations at e = 0.052, and at N — 256, 
we are able to establish coexistence of the two phases as 
in the monodisperse system Jig ]. Even there however, 
the N — > oo limit has remained controversial, and it is 
not yet firmly established whether the liquid-solid transi- 
tion in 2-d monodisperse disks is first-order or continuous 
(according to the KTNHY scenario, cf pH). Notice that 
the transition shows up very clearly both in the equation 
of state and in the compressibility, shown in the main 
figure Ha. 

In figure ||b, we show analogous data for e = 0.19. 
There, we see no indication of a thermodynamic phase 
transition in the whole range of densities studied. Sim- 
ilar results were obtained at e = 0.52 (remember that 
e < 0.58 cf 'text following eq.(ph); these data confirm and 
considerably extend our earlier calculations |Q. Within 
our flat probability distribution of the radii of particles, 
there exists a maximal polydispersity e — 0.58, which 
corresponds to the limit r m i n — ► in eq.(0). 

Our calculations thus confirm the existence of a ter- 
minal polydispersity, which was studied mostly in three 
dimensions, but also evoked in related two-dimensional 
work ill. 
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FIG. 3. Average bond-orientational order-parameter for 
P = 17.1. For polydispersities e < 0.08, the measure- 
ments converge towards a size-independent value, while we 
expect < <l/6 >^ for larger polydispersities. The proba- 
bility distributions of ^6 are unimodal, and give no indica- 
tion of a first-order phase transition. Notice the smoothness 
of the finite-size behavior, which is due to our choice of the 
fixed-probability increment model eq.(nl). 

The horizontal sweeps in the phase diagram of fig- 
ure [l] bracket the disordered-to-crystalline phase bound- 
ary, which must be essentially independent of e at large 
p. We further studied the thermodynamic system as 
a function of polydispersity at constant pressure, along 
the vertical paths indicated in figure [3. There, our data 
strongly favor a continuous transition: We did not find 
hysteresis in V(P), and detected no abrupt change of 
the volume distribution function in the transition region. 
We also computed the average bond-orientational order- 
parameter < v&g > (cf, e.g., p2] for definitions and a 
discussion). There also, the probability distribution of 
^6 has a single peak, and is perfectly reproduced, with- 
out any indications of hysteresis. We conclude that the 



transition is continuous. In figure || we show < \&6 > as 
a function of polydispersity e for different system sizes. 
A detailed, more rigorous study of the probable two-step 
melting process in this system is left for further work. 
Let us note that the numerical analysis of the standard 
Kosterlitz-Thouless transition has proven to be very sub- 
tle even in the prototypical XY- model p3| . 
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FIG. 4. Long-time diffusion constants as a function of par- 
ticle size TijL for two values of the polydispersity e. Left: 
e = 0.52, right: e = 0.19. In both cases, we witness a dramatic 
slow-down of the time-evolution as the density p ~ 0.8 is ap- 
proached. This density is much higher than the liquid-solid 
transition density of the monodisperse hard-disk system, but 
much smaller than the densities which we can reach (within 
thermal equilibrium) with the non-local algorithms. 

Finally, we determine the dynamic properties of the 
system, in the [NV] ensemble, with the local MC algo- 
rithm, and at values of e which ensure the absence of a 
liquid-solid transition. Using the protocol of ref. |7), we 
obtained the effective diffusion constants Di (cf figure |j), 
for e = 0.19 and e = 0.52. In both cases, the results are 
consistent with the scaling form f24|j: 



A(p) ~ (# - p) a , 



(4) 



where we extrapolate the diffusion constants for each par- 
ticle separately. The diffusion constants Di strongly de- 
pend on i, but the extrapolated transition densities pf 
do not, especially for small e. We take these extrapolated 
values as a definition of the glass transition density p . 
Only at large e (left part of figure) can some very small 
particles (r m i n /r max — 19) escape through slits left open 
in the system. This leads to some additional structure in 
the Di plot on the left half of figure 0, and to an increased 
value of the extrapolated pf for small i. 

In conclusion, our study of the polydisperse hard disk 
system has given the homogeneous phases shown in the 
schematic diagram figure [0. Inhomogeneous phases JO] 
seem to play no role at the parameters studied in this 
paper. A tendency towards phase separation should first 
show up in irregular finite-size behavior, which we have 
not encountered. We find it particularly significant that 
the glass transition depends so little on polydispersity, 
and takes place at densities (and pressures) much higher 



than the crystallization density in the monodisperse sys- 
tem. Using our specialized MC methods, we are able to 
probe the thermodynamics of the hard-disk system far 
within the glassy phase, where we find no indications of 
an accompanying thermodynamic transition. 
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